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ABSTRACT 

We use high-resolution N-body simulations to study the equilibrium density profiles 
of dark matter halos in hierarchically clustering universes. We find that all such 
profiles have the same shape, independent of halo mass, of initial density fluctuation 
spectrum, and of the values of the cosmological parameters. Spherically averaged 
equilibrium profiles are well fit over two decades in radius by a simple formula originally 
proposed to describe the structure of galaxy clusters in a cold dark matter universe. 
In any particular cosmology the two scale parameters of the fit, the halo mass and its 
characteristic density, are strongly correlated. Low-mass halos are significantly denser 
than more massive systems, a correlation which reflects the higher collapse redshift 
of small halos. The characteristic density of an equilibrium halo is proportional to 
the density of the universe at the time it was assembled. A suitable definition of 
this assembly time allows the same proportionality constant to be used for all the 
cosmologies that we have tested. We compare our results to previous work on halo 
density profiles and show that there is good agreement. We also provide a step-by-step 
analytic procedure, based on the Press-Schechter formalism, which allows accurate 
equilibrium profiles to be calculated as a function of mass in any hierarchical model. 
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1. Introduction 

It is twenty-five years since the discovery that galaxies are surrounded by extended massive 
halos of dark matter. A variety of observational probes - disk rotation curves, stellar kinematics, 
gas rings, motions of globular clusters, planetary nebulae and satellite galaxies, hot gaseous 
atmospheres, gravitational lensing effects - are now making it possible to map halo mass 
distributions in some detail. These distributions are intimately linked to the nature of the dark 
matter, to the way halos formed, and to the cosmological context of halo formation. 

Insight into these links came first from analytic studies. Building on the early work of Gunn & 
Gott (1972), similarity solutions were obtained by Fillmore Sz Goldreich (1984) and Bertschinger 
(1985) for the self-similar collapse of spherical perturbations in an Einstein-de Sitter universe. 
Such solutions necessarily resemble power laws in the virialized regions. Hoffman & Shaham (1985) 
and Hoffman (1988) extended this analysis by considering open universes, and by modeling as 
scale-free spherical perturbations the objects which form by hierarchical clustering from power-law 
initial density perturbation spectra (P(k) oc k n ). They argued that isothermal structure (p oc r~ 2 ) 
should be expected in an Einstein-de Sitter universe if n < —2, and that steeper profiles should be 
expected for larger n and in open universes. 

Despite the schematic nature of these arguments, their general predictions were verified as 
numerical data became available from iV-body simulations of hierarchical cosmologies. Power-law 
fits to halo density profiles in a variety of simulations all showed a clear steepening as n increases 
or the density of the universe decreases (Frenk et al. 1985, 1988; Quinn et al. 1986; Efstathiou et 
al. 1988; Zurek, Quinn & Salmon 1988; Warren et al. 1992; Crone, Evrard & Richstone 1994). 
An apparent exception was the work of West et al. (1987), who found that galaxy cluster density 
profiles show no clear dependence on n. 

Significant departures from power-law behaviour were first reported by Frenk et al. (1988), 
who noted that halo profiles in cold dark matter (CDM) simulations steepen progressively with 
increasing radius. Efstathiou et al. (1988) found similar departures - at odds with the analytic 
predictions - in their simulations of scale-free hierarchical clustering. They also noted that these 
departures were most obvious in their best resolved halos. Similar effects were noted by Dubinski 
& Carlberg (1991) in a high resolution simulation of a galaxy-sized CDM halo. These authors 
found their halo to be well described by a density profile with a gently changing logarithmic slope, 
specifically that proposed by Hernquist (1990). 

In earlier papers of this series we used high-resolution simulations to study the formation 
of CDM halos with masses spanning about four orders of magnitude, ranging from dwarf galaxy 
halos to those of rich galaxy clusters (Navarro, Frenk & White 1995, 1996). This work showed 
that the equilibrium density profiles of CDM halos of all masses can be accurately fit over two 
decades in radius by the simple formula, 

p( r ) _ $c m 
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where r s is a scale radius, S c is a characteristic (dimensionless) density, and p cr u = 3H 2 /8irG is 
the critical density for closure. This profile differs from the Hernquist (1990) model only in its 
asymptotic behaviour at r > r s (it tends to r~ 3 instead of r -4 ). Power-law fits over a restricted 
radial range have slopes which depend on the range fitted, steepening from —1 near the center to 
—3 at large r/r s . 

This similarity between CDM halos of widely differing mass is surprising in view of the strong 
dependence on power spectrum shape reported in earlier studies. The effective slope of the CDM 
power spectrum varies from n e ff ~ —2 on galactic scales to n c g » -1 on cluster scales, so one 
might have expected shallower profiles in galaxy halos than in clusters. In fact, the opposite is 
true; low-mass halos are denser, i.e. have higher values of S c , than high-mass halos. This property 
reflects the higher collapse redshift of the smaller systems. Power-law fits actually yield steeper 
slopes for less massive halos both when carried out at a fixed radius and when carried out at a 
fixed fraction of the virial radius (see Figures 3 and 4 of Navarro, Frenk & White 1996). 

We argue below that the apparent relationship between profile shape and initial power 
spectrum seen in earlier work results from systematic differences in the characteristic density of 
the halos chosen when comparing different models. This interpretation is reinforced by the recent 
work of Cole & Lacey (1996) and Tormen, Bouchet & White (1996), who find that the profiles 
of massive halos formed from power-law spectra are well described by eq. 1. They also confirm 
the strong correlation between <5 C and halo mass seen in our earlier work, and they find that, at 
a given mass, halos are denser when n is larger. Thus, the spectral index n seems to control the 
exact relationship between characteristic density and halo mass, rather than the effective slope of 
the density profile. 

It is clear from this discussion that a comprehensive study of this problem needs to consider 
the role of at least three factors: the halo mass, the power spectrum of initial density fluctuations, 
and the values of the cosmological parameters. In this paper we present the results of a large 
set of N-body simulations specifically designed to address these issues. We consider a variety of 
hierarchical clustering models, including CDM and power-law initial fluctuation spectra, as well 
as different values of the cosmological parameters Qo and A. In each cosmology we study halos 
spanning a large range in mass, carefully choosing numerical parameters so that all systems are 
simulated with comparable numerical resolution. 

The plan of the paper is as follows. Our numerical experiments are described in §2. In §3 we 
present our results and in §4 we discuss them in the context of earlier work. In §5 we summarize 
our main conclusions. An appendix lays out the formulae necessary to calculate analytically the 
density profile of an equilibrium halo of any mass in any hierarchical cosmology. 



2. The Numerical Experiments 
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2.1. The Cosmological Models 

We analyze the structure of dark matter halos in 8 different cosmologies. Five are Einstein-de 
Sitter (Qq = 1) models with various power spectra: the standard biased CDM spectrum (SCDM 
model: = 1, h = 0.5, <jg = 0.63) and four power-law or "scale-free" spectra with indices n = 0, 
—0.5, —1, and —1.5. Two further models also have power-law spectra (n = and —1) but in an 
open universe (Qq = 0.1). The last model we consider is a low-density CDM model with a flat 
geometry (CDMA: Qq = 0.25, A = 0.75, h = 0.75, as = 1.3). (Here and throughout this paper we 
express the cosmological constant A in units of 3H$, so that a low-density universe with a flat 
geometry has + A = 1. We also adopt the standard convention of writing the present Hubble 
constant as Ho = 100 h km s _1 Mpc -1 .) 

The normalization of the CDM models is specified by as, the rms mass fluctuation in spheres 
of radius 8 ft, -1 Mpc. Because of self-similarity, the normalization of the scale- free models is 
arbitrary. The evolutionary state of these models may be fully specified by a single parameter, 
the current value of the "nonlinear mass", M+(z). This mass scale is defined by requiring that the 
variance of the linear overdensity field at redshift z = 0, smoothed with a top-hat filter enclosing 
a mass M = M*, should equal the square of the critical density threshold for spherical collapse by 
redshift z: A(j(M*(z)) = b 2 crit {z, Qq,A). (See the Appendix for details on the computation of 6 cr u). 
This definition provides a "natural" way to scale the scale-free simulations to physical units and 
to compare different cosmological models. In scale- free models the mass scale defined by M*(z) is 
the only physical scale, and therefore the structure of halos can depend on mass only through the 
ratio M/M+. 

2.2. The Simulations 

Large cosmological N-body simulations are required to simulate the evolution of dark matter 
halos in their full cosmological context. However, such simulations are not generally well suited 
to explore a large range of halo masses. This is because systems of differing mass formed in 
a single simulation are resolved to differing degrees. More massive systems are better resolved 
because they contain more particles and because the gravitational softening is a smaller fraction 
of the virial radius. These systematic differences can introduce insidious numerical artifacts in 
the mass trends that we wish to investigate. We circumvent this problem by using the procedure 
outlined by Navarro, Frenk & White (1996). Halos are first identified in cosmological N-body 
simulations of large periodic boxes and then resimulated individually at higher resolution. During 
the resimulation the remainder of the original simulation is treated only to the accuracy needed 
to model tidal effects on the halo of interest. The advantage of this procedure is that numerical 
parameters can be tuned so that all halos are simulated with comparable resolution. Its main 
disadvantage is that only one halo is modeled per simulation, so that many simulations are needed 
to compile a representative halo sample. 
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2.2.1. The cosmological simulations 

The cosmological simulations were carried out using the P 3 M code of Efstathiou et al. 
(1985). The desired initial power spectrum was generated by using the Zel'dovich approximation 
to displace particles from a uniform initial load. The uniform load we used was either a "glass" 
configuration (White 1996) or a cubic grid. For simulations with power-law power spectra, the 
amplitude of the initial displacements was chosen by setting the power of the perturbed density 
field to the white-noise level at the Nyquist frequency of the particle grid. These simulations 
followed 10 6 particles on a 128 3 mesh and were stopped when the nonlinear mass, M*, corresponded 
to 1,000-2,000 particles. We identify this time with the present (z = 0). Since clustering evolves 
faster for more negative values of n, an expansion factor of 9.5 was sufficient for n = —1.5 whereas 
an expansion factor of 90 was required for n = (f2o = 1). Open models require even longer 
integrations, an expansion factor exceeding 150 for n = and f^o = 0.1. These simulations used 
the time-stepping and numerical scheme described in Efstathiou et al (1988). (Note, however, that 
the definition of in that paper differs from the one we use here.) 

We ran two P 3 M simulations for each of our CDM models. The SCDM runs followed 
64 3 particles in periodic boxes of 180 h~ x and 15 h^ 1 Mpc and were stopped when ag = 0.63 
(M* « 1.6 x 10 13 /i _1 Mo), the time which we identify with the present. The CDMA runs followed 
10 6 particles in boxes of 140 h' 1 and 46.67 h' 1 Mpc, until a 8 = 1.3 (M* « 4.1 x lO 13 /* -1 M ). 

2.2.2. The individual halo simulations 

Halos to be resimulated at higher resolution were selected randomly at z = from a list of 
clumps identified using a friends-of-friends group finder with linking length set to 10% of the mean 
interparticle separation. We chose masses in the range 0.1-lOM* for the power-law models and 
0.01-100M* for the CDM models. (The mass range is larger for the CDM models because in this 
case halos were chosen from two parent cosmological simulations with different box sizes.) Because 
we are interested the structure of equilibrium halos we were careful not to choose for analysis any 
halo which is far from virial equilibrium. In practice, we analyze each resimulated halo at the 
time between redshifts 0.05 and when it is closest to dynamical equilibrium, defined as the time 
when the ratio of kinetic to potential energy is closest to 0.5 for material within the virial radius. 
[Throughout this paper, we measure halo masses, M200, within a virial radius, r2oo, defined as the 
radius of a sphere of mean interior density 200 p cr i t . Halo circular speeds, V200 = (GM200 /^oo) 1 / 2 , 
are also measured at this radius unless otherwise specified. Numerical experiments show that for 
fi = 1 this radius approximately separates the virialized and infall regions (Cole & Lacey 1996). 
For convenience we continue to use these definitions when Qq ^ 1.] 

Once a halo is chosen for resimulation, the particles within its virial radius are traced back 
to the initial conditions, where a small box containing all of them is drawn. This box is filled 
with ~ 32 3 particles on a cubic grid which are then perturbed using the waves of the original 
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P 3 M simulation, together with extra high-frequency waves added to fill out the power spectrum 
between the Nyquist frequencies of the old and new particle grids. The regions beyond the 
"high-resolution" box are coarsely sampled with a few thousand particles of radially increasing 
mass in order to account for the large-scale tidal fields present in the original simulation. 

This procedure ensures the formation of a clump similar in all respects to the one selected in 
the P 3 M run, except for the improved numerical resolution. The size of the high-resolution box 
scales naturally with the total mass of each system and, as a result, all resimulated halos have 
about the same number of particles within the virial radius at z = 0, typically between 5, 000 
and 10, 000. The extensive tests presented in Navarro, Frenk & White (1996) indicate that this 
number of particles is adequate to resolve the structure of a halo over approximately two decades 
in radius. We therefore choose the gravitational softening, h g , to be 1% of the virial radius in all 
cases. (This is the scale-length of a spline softening; see Navarro & White 1993 for a definition.) 
The tree-code carries out simulations in physical, rather than comoving, coordinates, and uses 
individual timesteps for each particle. The minimum timestep depends on the maximum density 
resolved in each case, but was typically lO -5 //^ 1 . 

As discussed in Navarro, Frenk & White (1996), numerically convergent results require 
that the initial redshift of each run, Zi n u, should be high enough that all resolved scales in the 
initial box are still in the linear regime. In order to satisfy this condition, we chose z^nit so that 
the median initial displacement of particles in the high-resolution box was always less than the 
mean interparticle separation. Problems with this procedure may arise if Zi n a is so high that the 
gravitational softening (which is kept fixed in physical coordinates) becomes significantly larger 
that the mean initial interparticle separation. We found this to be a problem only for the smallest 
masses, M<M*, in the n = 0, Qq = 0.1 model. In this case the initial redshift prescribed by the 
median displacement condition is Zi n u > 700 and the gravitational softening is then a significant 
fraction of the initial box. This can affect the collapse of the earliest progenitors of these systems 
and so introduce spurious effects. We therefore limit our investigation to M>M+ in this particular 
cosmological model. Further tests of the effects of particle number, timestep size, and gravitational 
softening are given by Tormen et al. (1996). Their results confirm that the numerical parameters 
chosen here are adequate to give stable and accurate results. 



3. Results 

3.1. Time Evolution 

Figure 1 illustrates the time evolution of halos of different mass selected from the CIq = 1, 
n = — 1 series. Time runs from top to bottom and mass increases from left to right. The box size 
in each column is chosen so as to contain always approximately the same number of particles. The 
redshifts of each snapshot have been chosen so that the nonlinear mass M* increases by factors of 
4 from Z2 to Z\ and from z\ to zq = 0. This figure illustrates convincingly that low-mass halos 
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complete their assembly earlier than more massive systems. 



3.2. Density Profiles 

Figure 2 shows spherically averaged density profiles at z = for one of the least and one of 
the most massive halos for each set of cosmological parameters. These halos span almost four 
orders of magnitude in mass in the case of the CDM models, and about two orders of magnitude 
in mass in the power-law models. Radial units are kpc for CDM models (scale at the top), and are 
arbitrary in the power-law panels. Density is in units of lO lo M /kpc 3 in the CDM models and 
in arbitrary units in the others. Solid lines are fits to each halo profile using eq. 1. This simple 
formula provides a good fit to the structure of all halos over about two decades in radius, from the 
gravitational softening (indicated by arrows in Figure 2) to about the virial radius. The quality 
of the fit is essentially independent of halo mass or cosmological model and implies a remarkable 
uniformity in the equilibrium structure of dark matter halos in different hierarchical clustering 
models. 

The solid and dashed lines in Figure 3 show the profile fits of Figure 2 but with the radius 
scaled to the virial radius of each halo. This scaling removes the intrinsic dependence of size 
on mass (more massive halos are bigger) and allows a meaningful comparison between halos of 
different mass. From the definition of virial radius, the "concentration" of a halo, c = r2oo/ r s, and 
the characteristic density, 5 C , are linked by the relation 

200 c 3 

6c = ~ [ln(l + c) - c/(l + cjj ' (2) 

Thus at given halo mass (specified by M200), there is a single free parameter in eq. 1, which may 
be expressed either as the characteristic density, 5 C , or as the concentration parameter, c. This free 
parameter varies systematically with mass; Figure 3 shows that c and 5 C decrease with increasing 
halo mass. 

A universal density profile implies a universal circular velocity profile, V c (r) = (GM(r) /r) 1 / 2 . 
This is illustrated in Figure 4, where we plot V^-profiles for the same systems shown in Figure 2. 
As in Figure 3, radii are plotted in units of the virial radius; circular speeds have been normalized 
to the value at the virial radius, V2oo- The circular velocity curve implied by eq. 1 is 

fV c (r)\ 2 = 1 ln(l + cx) - (cx)/(l + ex) 

V V 200 ) x ln(l + c) - c/(l + c) ' U 

where x = r/r2oo is the radius in units of the virial radius. Circular velocities rise near the 
center, reach a maximum (Vrnax^) at x ma x ^ 2/c, and decline near the virial radius. More centrally 
concentrated halos (higher 5 C or higher c) are characterized by higher values of V max /V2oo. The 
dashed lines in Figure 4 show plots of eq. 3 with parameter values derived from the fits to the 
density profiles of Figure 2. The dotted lines are fits using a Hernquist (1990) model constrained 
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to match the location of the maximum of the V^-curve. The two fits are indistinguishable near 
the center, but the Hernquist model underestimates V c near the virial radius. This disagreement 
becomes more pronounced in lower mass systems, for which 5 C and V max /V2oo are larger. 



3.3. The mass dependence of halo structure 

The mass-density dependence pointed out above is further illustrated in Figure 5, where we 
plot 5 C versus mass (expressed in units of M*) for all the systems in each series. An equivalent 
plot, illustrating the mass dependence of the concentration, c, is shown in Figure 6. (The panel 
on the upper-left, corresponding to the SCDM model, is equivalent to Figure 7 of Navarro, Frenk 
& White 1996). The characteristic density of a halo increases towards lower masses in all the 
cosmological models considered. This result supports the idea that the M2qq-5 c relation is a direct 
result of the higher redshift of collapse of less massive systems, and suggests a simple model to 
describe the mass-density relation. This model assigns to each halo of mass M (identified at z = 0) 
a collapse redshift, z co u(M,f), defined as the time at which half the mass of the halo was first 
contained in progenitors more massive than some fraction / of the final mass. With this definition, 
Zcoii can be computed simply using the Press-Schechter formalism (e.g. Lacey & Cole 1993), 

erfc( , S crit(Zcoll) ~ S°crit ) = ^ (4) 



^(AgtfMj-AgCM))/ 2 ' 



where Ag(M) is the linear variance of the power spectrum at z = smoothed with a top-hat filter 
of mass M, 5 cr it(z) is the density threshold for spherical collapse by redshift z, and 5® rit = 5 cr u(0). 
[This definition can be extended to halos identified at any redshift z by replacing 5® rit by S cr it(zo) 
in eq. 4.] Assuming the characteristic density of a halo to be proportional to the density of the 
universe at the corresponding z co u then implies 

5 c (M\f) = Cn (l + Zcoii(MJ)f (5) 

where C is a proportionality constant which might, in principle, depend on / and on the power 
spectrum. 

We will see below that / <C 1 is needed for this argument to give a good fit to our simulation 
data. In this limit Aq(/M) > Aq(M) and eq. 4 reduces to 

ScHt(zcoii) = S° cnt + C'Ao(fM), (6) 

where C « 0.7. For / < 1, 5 cr it(z co u) 3> 5® rit for all masses in the range of interest, so that 
Scrit(zcoll) « A (/M). Since M±(z coU ) is defined by A (M^(z co a)) = S cri t(z co ii), this equation 
implies that the characteristic density of a halo is proportional to the mean density of the universe 
at the time when M* « fM, ie. when the characteristic non-linear mass is a fixed small fraction 
of the final halo mass. For scale-free models this implies <5 C oc M - (™ +3 )/ 2 , the same scaling that 
links M*(z) and the mean cosmic density at redshift z. 
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Figure 5 shows the correlations predicted from eq. 5 for three values of the parameter /: 0.5, 
0.1, and 0.01. The value of the proportionality constant, C(/), is chosen in each case in order to 
match the results of the Einstein- de Sitter simulations for M = M*. These values are given in 
Table 1. The same values of C{f) are used to plot the curves in the panels corresponding to the 
low-density models. Some interesting results emerge from inspection of Figure 5 and Table 1. 

(i) The agreement between the mass-density dependence predicted by eq. 5 and the results of 
the Einstein-de Sitter simulations improves for smaller values of /. This is also true for the 
low-density models. Once C(f) is fixed by matching the results of the Einstein-de Sitter 
models, the same value of C(f) provides a good match to the low density models only if /<0.01. 
Interestingly, for / = 0.01 approximately the same value of the proportionality constant, 

C 3 x 10 3 , seems to fit all our simulations. 

(ii) The characteristic density of halos decreases systematically for more negative values of 
the spectral index n. At M = M±, SCDM halos are the least dense in our f2o = 1 series, less 
concentrated still than those corresponding to n = —1.5. This is consistent with the general 
trend because, according to eqs. 4 and 5, the characteristic density of a halo of mass M± is 
controlled by the shape of the power spectrum on scales ~ fM±. This is about ~ 10 11 M Q for 
/ f» 0.01 and the effective slope of the CDM spectrum on this mass scale is n e g ~ —2. 

(iii) For the power-law models with n = and —1 the characteristic density at a given M/M* 
increases as Qq decreases. Such a trend is plausible since we expect the collapse redshift of 
halos of a given mass to increase as Qq decreases. On the other hand, halos formed in the 
low-density CDMA universe are actually less dense than those formed in the standard biased 
CDM model because S c depends not only on collapse redshift but also on f^o (see eq. 5). 
Although reducing f^o increases the collapse redshift, the increase in 5 C from the (1 + z co u) 3 
factor can be outweighed by the change in Qq. In the CDMA model the two effects can combine 
to give a reduction in 5 C as Qq decreases. (We remind the reader that 5 C is defined relative to 
the critical density rather than the mean density.) 

(iv) Each halo has a characteristic maximum circular speed, V max (see eq. 3), which is strongly 
correlated with its mass. This is shown in Figure 7, where we also plot least squares fits of the 
form M200 Vmax to the data in each series. Consistent with the trends shown in Figures 
5 and 6, the correlation steepens as n increases; we find a ~ 3.2 - 3.3 for CDM models and 
a > 5 for n = 0. Note also that the correlations are extremely tight; the rms scatter in 
logM2oo is less than ~ 0.1 in all cases. This is in part a consequence of the generally good 
fit of the density profile of eq. 1. The ratio V max /V200 increases only logarithmically with the 
central concentration of the halo, and changes only by about a factor of two as 5 C varies by 
four orders of magnitude between 10 3 and 10 . As a result, the M2oo-V max relation does not 
deviate much from the M200-V200 relation which, by definition, has zero scatter. This has 
important consequences for the expected tightness of empirical correlations between mass and 
characteristic velocity, such as the Tully-Fisher relation. We intend to return to this issue in 
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future work. 

The results in Figures 5-7 support our conclusion that the characteristic density of a halo 
is controlled mainly by the mean matter density of the universe at a suitably defined time of 
collapse. One important test of this interpretation is to measure z co u directly in the simulations 
and to compare the result to eq. 5. To do this we identify clumps at every output time using our 
friends-of-friends group-finder with linking length set to 20% of the current mean interparticle 
separation. We then trace the particles in the most massive clump identified at z = (which 
typically has a mean overdensity of ~ 200) and, at each redshift, add up the total mass in clumps 
which contain any of these particles and which are individually more massive than 10% of the final 
mass. We identify z co u with the redshift at which this mass first exceeds half of the final mass. 
This is roughly equivalent to the analytic procedure outlined in eq. 4 for / = 0.1. (We decided to 
use / = 0.1 rather than / = 0.01 because the smaller value results in very high collapse redshifts, 
often before the first output in the simulation.) The main difference is that some of the mass from 
the high redshift progenitors ends up outside the virial radius at z = 0. This causes a slight bias 
of the measured collapse redshifts towards higher values than the Press-Schechter predictions. 

The correlation between d c /Qo and (1 + z co u) 3 obtained using this procedure is shown in 
Figure 8. All halos in Einstein-de Sitter models are shown with filled circles; those in open 
universes with open circles; and those in the CDMA model with starred symbols. The solid line is 
the relation predicted by eq. 5 for C = 5 x 10 3 . This is clearly an excellent approximation to the 
results of the numerical simulations, and confirms that the mean matter density of the universe 
at the time of collapse is the main factor determining the characteristic density of a halo. Note 
that the value of the proportionality constant is slightly lower than the values given in Table 
1 for / = 0.1. This difference compensates for the slight bias towards higher collapse redshifts 
introduced by our numerical procedure. 

A summary of our results is presented in Figure 9. The panels on the left compile the fits 
(for / = 0.01) to the mass dependence of 5 C and c in the Einstein-de Sitter models. The panels 
on the right compare the Einstein-de Sitter results with those for low-density models. As noted 
above, the typical density of M*-halos increases with n. However, the difference between models 
becomes less pronounced at higher masses, and is almost negligible at M > 10M*. Halos of a given 
M/M+ in low density universes can have either lower or higher characteristic densities than their 
Einstein-de Sitter counterparts, depending on the competing effects of the collapse redshift and 
the value of Note that all the curves in Figure 9 use the same value, / = 0.01, and essentially 
the same value of the proportionality constant C ~ 3 x 10 3 (see eq. 5). Thus, once calibrated 
for an Einstein-de Sitter model, it is possible to apply eqs. 4 and 5 to predict the characteristic 
density of halos formed in other hierarchically clustering models. In the Appendix we provide a 
detailed description of how to compute 5 C (M) numerically for a variety of cosmologies. 



- 11 - 



3.4. The scatter in the correlations 

We examine now the origin of the scatter in the correlations presented above. In particular, 
we explore whether at fixed mass the dispersion in the measured values of S c is due to variations 
either in the collapse redshift or in the global angular momentum of the system. As shown by 
Figures 10 and 11, the bulk of the scatter in 5 C at a given M/M+ can be attributed to small 
differences in the redshift of collapse. Figure 10 shows that the (5 c -deviations from the solid-line fits 
in Figure 5 (i.e. those for / = 0.01) correlate strongly with deviations in the redshift of collapse. 
(The latter are measured from least-square fits to the M200 vs (1 + z co u) correlations measured 
directly from the simulations.) Furthermore, the two residuals seem to correlate just as expected 
from eq. 5, i.e. Alog<5 c = 3Alog(l + z co u). (We note that the magnitude of this scatter may not 
be fully representative of the dispersion in halo properties corresponding to each cosmological 
model, since the sample has been selected so as to minimize departures from equilibrium.) This 
relation is indicated by the solid line and is seen to reproduce very well the trend observed in 
Figure 10. Figure 11 shows the same <5 c -residuals of Figure 10, but now as a function of the scatter 
in the mass-rotation parameter (M200 vs A) correlation. (The rotation parameter is defined by 
A = J\E\ 1 / 2 /GM^I 2 , where J is the angular momentum and E is the binding energy of the halo. 
The median A in our series is ~ 0.04, in good agreement with previous studies.) We find no 
discernible correlation between M/M* and A or between the 5 C - and A-residuals (Figure 11). This 
provides further evidence supporting our contention that the redshift of collapse is the primary 
factor determining the characteristic density of a halo. 

4. Comparison with previous work 

The main conclusion of our study, that the shape of halo density profiles is independent of 
cosmological context, appears to contradict previous work on this subject. As discussed in §1, 
a strong dependence of the slope of the density profile on the spectral index n and the density 
parameter has been established by a number of analytic and numerical studies. We now show 
that our results actually include and extend those of previous workers, and we offer an attractive 
explanation for some discrepancies found in the literature. 

We first address the claim by Quinn et al. (1986), Efstathiou et al (1988), Zurek et al. (1988), 
and Warren et al. (1992) that halo density profiles steepen with increasing values of the spectral 
index n in an Einstein-de Sitter universe. This claim is based on the discovery that circular 
velocity profiles of galaxy-sized halos are relatively flat for n ~ — 2 (or for SCDM; Frenk et al. 
1985) but decline progressively faster at large r for larger n. Figure 12 shows that we find the 
same trend if we analyze our data in the same fashion as these authors. In this figure we plot, in 
linear units, the VVprofile of an M + halo for different values of n. Each curve has been computed 
using eq. 3 and the values of c obtained from the fits presented in Figure 9. Since all these halos 
have the same mass, they also have the same virial radius and circular speed, r2oo and V^ooj 
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respectively, which we have used as normalizing factors. The linear units in Figure 12 obscure the 
fact that the form of the curves is the same in all cases, and encourage one to conclude, as did 
previous authors, that halo rotation curves steepen with increasing n. In fact this trend is present 
because M* halos collapse earlier, and so are denser, for larger values of n. 

We note that the trend with n in Figure 12 is sensitive to the choice of halo mass. Had we 
used halos with M < M* the trend would have been stronger. On the other hand, very massive 
halos (M >10 M*) have similar characteristic densities, irrespective of n (see Figure 9). Thus, 
had we plotted very massive halos in Figure 12 we would have concluded that the density profile 
depends only weakly on n (at least for VLq = 1). This is reminiscent of the claim by West et al. 
(1987) that the structure of galaxy cluster halos is independent of n in Einstein-de Sitter universes. 
Our analysis suggests an explanation for this apparently discrepant result. By considering only 
the most massive systems in their simulations, West et al. focused on a mass range where the 
dependence of 5 C on n is minimal, and so concluded (correctly) that cluster profiles depend very 
weakly on initial power spectrum. 

A similar explanation accounts for the weak dependence of halo density profile on spectral 
index n and on VLq reported by Crone et al. (1994). These authors chose to combine the 35 most 
massive clumps in each of their cosmological simulations in order to produce an "average" density 
profile for each value of n and Qq. They then fitted power laws of the form p(r) oc r 7 to these 
profiles in the radial range corresponding to density contrasts between 100 and 3000. For f2o = 1, 
this procedure yielded 7 values that decrease from —2.2 to —2.5 as n increases from —2 to (see 
the curve labeled "EdS" in their Figure 4). Our results show this weak trend to be a consequence 
of the averaging procedure they adopted. The shape of the halo mass function depends strongly 
on n and is increasingly peaked around M « M* for larger values of n. Thus, combining the 35 
most massive clumps in a simulation results in different "effective" masses for different values of 
n, with a bias towards larger values of M/M+ for more negative values of n. Applying the same 
selection procedure as Crone et al. to our own cosmological simulations we find that the median 
mass of these halo ensembles increases from ~ M* for n = to ~ 6M + for n = —1.5. As illustrated 
in Figure 13, fitting power-laws to the density profiles of these "average" halos results in a weak 
steepening with n similar to the trend reported by Crone et al. The slopes of 7 ~ —3 which they 
found for low density models are also easily understood, since in these cases they fit a radial range 
which extends well outside our nominal virial radius r2oo- 

Our results are also in agreement with recent work by Cole & Lacey (1996) and Tormen et 
al. (1996), who analyzed the structure of massive halos (M > M*) in scale-free universes. The 
general correlations with mass which they report (more massive halos tend to be less centrally 
concentrated than less massive halos) agree well with the results we present in §3. At a given halo 
mass, however, Cole & Lacey's halos are significantly less concentrated than ours. For example, 
for M + halos they find c ~ 12 and ~ 17 for n = — 1 and 0, respectively (f2o = 1), compared with 
c ~ 17 and ~ 30 in our simulations (see connected symbols in the lower-right panel of Figure 9). 
On the other hand, the results of Tormen et al. (1996) are in very good agreement with ours; 
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they find c ~ 10 for a 10M* halo formed in an n = — 1, J7o = 1 universe, which compares well with 
c ~ 9 that we find for a similar system. 

As discussed in detail by Tormen et al, the discrepancy between our results and those of Cole 
& Lacey is likely to be due in part to the poorer numerical resolution of their simulations. Indeed, 
their simulations used gravitational softenings that are 2-3 times larger than ours and timesteps 
which are also significantly longer than the typical values in our simulations. Both of these 
effects can artificially lower the central concentration of a halo. A concurrent factor may be the 
averaging procedure adopted by Cole & Lacey. These authors constructed average density profiles 
by co-adding all halos of similar mass identified in their cosmological simulations, regardless of 
their dynamical state. This sample contains a number of unrelaxed systems and ongoing mergers 
where substructure and double centers can bias the results of the fitting procedure towards lower 
concentrations. 

We can test directly for the effects of these various factors by applying the averaging procedure 
of Cole & Lacey to halos identified in our own P 3 M cosmological simulations (see §2.2.1) and 
comparing with the results of the individual halo runs. We restrict this comparison to the most 
massive halos in each run, M>10M^, since these systems have a large number of particles 
and comparable numerical resolution to that of Cole & Lacey. The comparison confirms that 
the central concentration of halos can be significantly underestimated as a result of the factors 
mentioned above. The magnitude of the effect is sensitive to n and f^oi the concentration c can be 
underestimated by up to a factor of ~ 3 for n = and Qq = 0.1, but by less than ~ 1.5 for n = — 1 
and Qq = I. We conclude that the disagreement between our results and those of Cole & Lacey is 
likely to be the result of the combined effect of their halo selection and averaging procedures and 
their poorer numerical resolution. 

5. Discussion 

Our results suggest that equilibrium dark halos formed through dissipationless hierarchical 
clustering have density profiles with a universal shape which does not depend on their mass, 
on the power spectrum of initial fluctuations, or on the cosmological parameters and A. It 
appears that mergers and collisions act during halo formation as a "relaxation" mechanism to 
produce an equilibrium which is largely independent of initial conditions. This mechanism must 
operate rapidly since the similarity between profiles extends to the virial radius. These properties 
are characteristic of the "violent relaxation" process proposed by Lynden-Bell (1967) to explain 
the regularities observed in the structure of elliptical galaxies, and it is interesting that our 
universal profile is similar to the Hernquist profile which gives a good description of elliptical 
galaxy photometry. The two differ significantly only at large radii, perhaps because ellipticals 
are relatively isolated systems whereas dark halos are not. Syer & White (1996) suggest that a 
universal profile may be understood as a fixed point for a process of repeated mergers between 
unequal objects. Their analysis of this process predicts a dependence on initial power spectrum 
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which seems stronger, however, than that seen in our numerical data. 

Our simulations suggest that the density profile of an isolated equilibrium halo can be 
specified quite accurately by giving two parameters; halo mass and halo characteristic density. 
Furthermore, in any particular hierarchical model these parameters are related in such a way that 
the characteristic density is proportional to the mean cosmic density at the time when the mass of 
typical nonlinear objects was some fixed small fraction of the halo mass. The characteristic density 
thus reflects the density of the universe at the collapse time of the objects which merge to form 
the halo core. With this interpretation we are able to fit the mass-density relation for equilibrium 
halos in all the cosmological models we have considered. In addition, it is possible to calculate the 
mass-density relation in any other hierarchical model (see Appendix). It is difficult to imagine 
a simpler situation - halos of all masses in all hierarchical cosmologies look the same, and their 
characteristic densities are just proportional to the cosmic density at the time they "formed". 

Our results extend the radial range over which dark halo structure can reliably be determined 
to more than two decades. The central regions of our models have densities of order W 6 p cr i t , 
comparable to those in the luminous parts of galaxies and in the central cores of galaxy clusters. 
As a result a variety of direct observational tests of our predictions are available. In Navarro 
et al. (1996) we discussed several of these in the context of the standard CDM cosmogony - 
rotation curves of giant and dwarf galaxies, satellite galaxy dynamics, hot gaseous atmospheres 
around galaxies and in clusters, strong and weak gravitational lensing - all provide interesting 
constraints. We will not pursue these issues here, however, because they would take us too far 
from our primary goal, namely the presentation of a simple and apparently general theoretical 
result: hierarchical clustering leads to a universal halo density profile just as it leads to universal 
distributions of halo axial ratios and halo spins; none of these properties depends strongly on 
power spectrum, on Q, or on A [see Cole & Lacey (1996) and references therein]. Comparison 
of the predicted halo structure with observation should, of course, provide strong constraints on 
the parameters which define particular hierarchical cosmogonies, and perhaps on the hierarchical 
clustering paradigm itself. We expect to come back to these issues in future work. 

We are grateful for the hospitality of the Institute for Theoretical Physics of the University 
of California at Santa Barbara, where some of the work presented here was carried out. JFN is 
also grateful for the hospitality of the Max Planck Institut fur Astrophysik in Garching, where 
this project was started. In addition he would like to acknowledge useful discussions with Cedric 
Lacey. Special thanks are due to Shaun Cole for making the data shown in Figure 9 available in 
electronic form. This work was supported in part by the PPARC and by the NSF under grant 
No. PHY94-07194 to the Institute for Theoretical Physics of the University of California at Santa 
Barbara. 
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A. Step-by-step calculation of the density profile of a dark matter halo 

In this Appendix we describe in detail the calculation of the parameters that specify the 
density profile of a dark halo of mass M. This calculation is applicable to Einstein-de Sitter 
(Qq = 1, A = 0), open (Qq < 1, A = 0), and flat (CIq + A = 1) universes. We provide approximate 
fitting formulae that are valid for power-law and CDM initial density fluctuation spectra. 



A halo of mass M identified at z = zq can be characterized by its virial radius, 

' m y/y n 
h-iM e ) \n(z ) 



-\/r \ 1/3 / Q \ —1/3 

r 200 = 1.63 xl(n 2 (— — - — °- (l + zd^/i^kpc, (Al) 



or by its circular velocity, 



gmV 2 ( r 200 \( n \V2. 



The density profile of this system is fully specified by its characteristic density <5 C and is given by 
(see eq.l) 

P(r) = £3l (1 + ^o) 3 ^ *° , 2 , (A3) 

8lT& ii{Zo) CX{1 + cx) z 

where x = r/r2oo and c is the concentration parameter, a function of S c given in eq. 2. The 
corresponding circular velocity profile, V c (r), is given by 

V c {r)\ 2 _ 1 ln(l + cx) - (cx)/(l + cx) 



V 2 oo / x ln(l + c) - c/(l + c 

using the concentration c as a parameter. 

The characteristic density is determined by the collapse redshift z co u, which is given by (see 
eq. 4) 

) S° crit (n(z coll ),A) D(z Q ,n ,A) n , 0.477 



1 + ^-TTT P(^IUM) - Ag(M)) (A5) 



^crit(^o) 6° rit (n(z ),A) D(z coU ,Q ,A) 6 crit (z ) 
Here D(z,r2o,A) is an ^-dependent linear growth factor that can be written as 



D(z,n ,A) = { 



'l/(l + z) if O = 1 and A = 0, 

F x (w) I F x (w ) if n < 1 and A = 0, (A6) 

. F2(y)F 3 (y)/F 2 (y )F 3 (yo) if n + A = 1, 



where we have used the following auxiliary definitions, 



w = +--l, (A7) 
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and 



= 1 + - + 3(1 \f 2 In [(1 + „)V2 _ n V2 L (ilg) 

y = (2u; ) 1/3 , {Aid) 

v = T+~z' M 

(u 3 + 2) 1 ! 2 

F 2 (u) = { + /2 j , (A12) 

/■«/ y' \ 3 / 2 



V 3 + 2 

A good numerical approximation to the O-dependence of the critical threshold for spherical 
collapse is given by 



<&i t (n,A) 



' 0.15 (12vr) 2 / 3 if f} = 1 and A = 0, 

0.15 (12vr) 2 / 3 ft - 0185 if fi < 1 and A = 0, (AU) 
k 0.15(127r) 2 /3QO.O055 ifn + A = l, 



which can be used to compute 5 cr it(zo) = 6^ rit (Q(zo))/D(zo, Qo, A). Finally, eq. A5 requires 
Ag(M), the variance of the power spectrum on mass scale M, extrapolated linearly to z = 0. In 
the case of a power-law spectrum of initial density fluctuations, P(k) oc k n , this is simply 

/ M \ -("+ 3 )/ 6 
Al(M) = S° crit — - - , (A15) 



M*(z = 0) 

where we have normalized the spectrum by M±{z = 0), the present nonlinear mass. A CDM 
spectrum is usually normalized by a&, the rms mass fluctuations within a sphere of radius 8/i _1 
Mpc, and its variance can be approximated by, 

A (M) = a 8 F 4 {M 8 )/F 4 (M h ), (A16) 

where we have used the following definitions, 

M s = 6.005 x 10 14 (/if} ) 3 , (-417) 



and 

F A {u) = A lU 0S7 [l + {A 2 u~ G - 1 + A 3 u-° m ) p ] l/p , (A19) 

with Ai = 8.6594 x 10" 12 , A 2 = 3.5, A 3 = 1.628 x 10 9 , and p = 0.255. 

Eqs. A6-A19 can be used to solve eq. A5 and find the collapse redshift, z co u, corresponding to 
a halo of mass M. As noted when discussing Fig. 5, we recommend using / = 0.01 when solving 
eq. A5, since this value seems to reproduce well the results of all our numerical experiments. 
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Once Zcoii(M) has been found, the characteristic density of the halo (expressed in units of the 
critical density at z = zq) can be computed from (see eq.5 and Table 1), 

S C (M, zo) ~ 3 x 10 3 n(z ) ( L±^*L) 3 , (A20) 

V t + zo J 

where we have assumed / = 0.01. A FORTRAN subroutine that implements the procedure 
described here and returns S C (M, zq) for all the cosmologies we discuss is available from the authors 
upon request. 
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P(k) 


n 


A 
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C(f) 


CDM 


1.0 


0.0 


0.5 


1.75 


x 10 4 




1.0 


0.0 


0.1 


7.44 


x 10 3 
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0.01 


3.41 


x 10 3 


CDMA 


0.25 


0.75 


0.5 
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x 10 4 




0.25 
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x 10 3 
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0. 
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X 


10 5 
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0.01 
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X 


10 3 



Table 1: Parameters in eq. 5 used to plot the fits in Figure 5. 
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M < M, M ~ M, M > M, 

Fig. 1. — Particle plots illustrating the time evolution of halos of different mass in an = 1, 
n = — 1 cosmology. Box sizes of each column are chosen so as to include approximately the same 
number of particles. At zq = the box size corresponds to about 6 x r2oo- Time runs from top to 
bottom. Each snapshot is chosen so that M* increases by a factor of 4 between each row. Low mass 
halos assemble earlier than their more massive counterparts. This is true for every cosmological 
scenario in our series. 
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Fig. 2. — Density profiles of one of the most and one of the least massive halos in each series. In 
each panel the low-mass system is represented by the leftmost curve. In the SCDM and CDMA 
models radii are given in kpc (scale at the top) and densities are in units of lO lo M /kpc 3 . In all 
other panels units are arbitrary. The density parameter, Qq, and the value of the spectral index, n 
is given in each panel. Solid lines are fits to the density profiles using eq. (1). The arrows indicate 
the value of the gravitational softening. The virial radius of each system is in all cases two orders 
of magnitude larger than the gravitational softening. 
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Fig. 3. — The fits to the density profiles of Figure 2, scaled to the virial radius, r^oo, of each 
system and to the critical density of the universe at z = 0. Solid and dashed lines correspond to the 
low- and high-mass systems, respectively. Note that low-mass systems are denser than high-mass 
systems near the center, indicating that the characteristic density of a halo increases as the halo 
mass decreases. 
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Fig. 4. — The circular velocity profiles of the halos shown in Figure 2. Radii are in units of the 
virial radius and circular speeds are normalized to the value at the virial radius. The thin solid line 
shows the data from the simulations. All curves have the same shape: they rise near the center 
until they reach a maximum and then decline at the outer edge. Low mass systems have higher 
maximum circular velocities in these scaled units because of their higher central concentrations. 
Dashed lines are fits using eq.(3). The dotted lines are the fit to the low-mass halo in each panel 
using a Hernquist profile. Note that this model fits rather well the inner regions of the halos, but 
underestimates the circular velocity near the virial radius. 
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Fig. 5. — The correlation between the mass of a halo and its characteristic density. Masses are 
given in units of the nonlinear mass scale M* (see text for a definition). Densities are relative to 
the critical value. Three curves are shown in each panel for different values of the parameter / (see 
eq. 5). The fits are normalized to intersect at M200 = M* in the case Q = 1. This normalization 
is then used for the low-density models (Q,q < 1). Note that for / = 0.01 this procedure results in 
good fits to the results of the simulations in all cases. 
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Fig. 6.— 



As Figure 5, but for the concentration parameter c. 
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Fig. 7. — The mass dependence of the maximum circular velocity of a halo. Velocity units are 
arbitrary in the power-law panels and km s _1 in the CDM models (scale at top). Mass is in units 
of 1O 1O M0 for CDM halos and in units of M + in the other panels. Power law fits of the form 
M oc V^ax are shown. The value of a and the rms scatter in the mass about the fit are indicated 
in each panel. Note that the M-V max dependence steepens for larger values of the spectral index 
n. The effect of the cosmological parameters on a seems to be rather small. 
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Fig. 8. — The characteristic density of all halos in our series as a function of the redshift at which 
half of the final mass is in collapsed progenitors more massive than 10% of the final mass. Solid 
circles correspond to all our runs with 0,q = 1, open circles to our runs with Qq = 0.1 and A = 0, 
and starred symbols to the CDMA runs (f^o = 0.25, A = 0.75). The solid line shows the "natural" 
scaling, 5 C oc f&o(l + z) 3 expected if the characteristic density of a halo is directly proportional to 
the mean matter density of the universe at the time of collapse. 
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Fig. 9. — A summary of the mass dependence of 5 C and c in our different cosmological models. 
Mass is given in units of the nonlinear mass scale M*. Curves are labeled in the upper plots and 
the same line types are used in the bottom plots. The symbols in the lower-right panel show the 
correlation between c and mass found by Cole & Lacey (1996), with open squares for n = — 1, 
fio = 1 and solid squares for n = 0, f2o = 1. These results should be compared to the lower dashed 
and lower dot-dashed curves in the same panel, respectively. Because of their poorer numerical 
resolution, Cole & Lacey's halos are significantly less concentrated than the ones in our study. 
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-0.4 -0.2 0.2 0.4 

Alog(l + z coU ) 

Fig. 10. — The <5 c -deviation from the fits shown in Figure 5 (for / = 0.01) plotted versus the 
deviation from the mean relation between the collapse redshift, l + z co u (Figure 8), and the mass of 
a system. Note that systems assembled earlier (later) than the average tend to have characteristic 
densities above (below) the mean. The correlation between the 6 C and 1 + z co u residuals follows 
the "natural" scaling, Alog<5 c = 3Alog(l + z co u), as shown by the solid line. Only the results 
corresponding to Qq = 1 are shown. 
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Fig. 11. — The <5 c -residuals in Figure 5 relative to the solid curve fit (/ = 0.01) plotted versus the 
residuals in the M200-X correlation. This figure shows that the scatter in 8 C at a given mass cannot 
be attributed to the effects of rotation. All halos with = 1 are included in this plot. 
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Fig. 12. — Circular velocity profiles of M± halos formed in simulations with different power spectra 
in an Einstein-de Sitter universe. The curves were computed using eq. 3 and the values of the 
concentration c obtained from the fits in Figure 9. Radii are in units of the virial radius and 
circular speeds in units of the circular velocity at the virial radius. Note that, because of the use 
of linear units, the fact that all curves have the same shape (eq. 3) is not immediately apparent. 
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Fig. 13. — A comparison between our density profiles and the power-law fits of Crone et al. (1994). 
The region fitted by these authors, corresponding to densities between 100 and 3000, is shown by 
the dotted lines. Over this region our results (solid lines) and the power-law parameterization 
adopted by Crone et al. (dashed lines) are consistent. The profiles shown by the solid curves 
correspond to M = 6M* (n = -1.5), M = 4M* (n = -1), M = 2M* (n = -0.5), and M = M* 
(n = 0), as discussed in the text. 



